#ifndef _MVN_H_
#define _MVN_H_

#include <vector>
#include <boost/math/special_functions/erf.hpp>

using namespace std;

template<typename T>
class Matrix
{
public:
	//Constructor using the matrix size
	Matrix(size_t size):m_size(size),m_mat(m_size)
	{
		for(typename vector<vector<T> >::size_type i=0; i<m_size; ++i)
		{
			m_mat[i] = vector<T>(m_size);
		}
	}

	//Constructor using matrix values
	Matrix(T* narray, size_t t): m_size(t),m_mat(m_size)
	{
		for(size_t i=0; i<t; i++){
			for(size_t j=0; j<t; j++){
				m_mat[i].push_back(narray[i*t+j]);
			}
		}
	} 
    
	//Undefined behavior if index underflow or overflow.
	const vector<T>& operator [] (size_t row) const
	{
		return m_mat[row];
	}
	vector<T>& operator [] (size_t row)
	{
		return m_mat[row];
	}

    size_t getSize() const
	{
		return m_size;
	}
	
private:
	size_t					m_size;
	vector<vector<T> >		m_mat;
};


//Auxiliary functions

//Calculate the Cholesky decomposition of the matrix 'mat'
Matrix<double> CholeskyFactor(const Matrix<double> &mat);

//Error function which is used in the function 'Phy'
inline double ErrFunc(double x) 
{
	return boost::math::erf(x);
}

//Inverse error function
double InvErrFunc(double);

inline double Phy(double x)
{
	return 0.5 + 0.5 * ErrFunc(x/sqrt(2.0));
}

inline double InvPhy(double x)
{
	return sqrt(2.0)*InvErrFunc(2*x-1);
}

inline int Sgn(double x)
{
	return x>0?1:(x<0?-1:0);
}

#endif